##########################
# Script to reproduce the figures and results in Busemeyer & Tober: 
# "Dealing with Technological Change: Social Policy Preferences and 
# Institutional Context"
# Published in Comparative Political Studies
# October 2022
##########################

##########################
# Load packages
##########################

# MCMCglmm (required for multivariate mixed-effects models)
if (!require("MCMCglmm")) install.packages("MCMCglmm")
library(MCMCglmm)
# arm (required for standardizing data)
if (!require("arm")) install.packages("arm")
library(arm)
# plyr (required for recoding via mapvalues)
if (!require("plyr")) install.packages("plyr")
library(plyr)
# tidyr (required for graphs)
if (!require("tidyr")) install.packages("tidyr")
library(tidyr)
# dplyr (required for graphs)
if (!require("dplyr")) install.packages("dplyr")
library(dplyr)
# broom.mixed (required for graphs)
if (!require("broom.mixed")) install.packages("broom.mixed")
library(broom.mixed)
# dotwhisker (required for graphs)
if (!require("dotwhisker")) install.packages("dotwhisker")
library(dotwhisker)
# ggplot2 (required for graphs)
if (!require("ggplot2")) install.packages("ggplot2")
library(ggplot2); theme_set(theme_bw())
# haven (required for dealing with labels)
if (!require("haven")) install.packages("haven")
library(haven)

##########################
# Load data
##########################

load("dat.RData")

##########################
# Figure 1
##########################

dat.des <- as.data.frame(dat[c("train","retrain","redis","ubi","subs",
                               "country","saut")])
policies_des <- dat.des %>%
  group_by(country) %>%
  summarize_at(vars(c("train","retrain","redis","ubi","subs")),
               list(name=mean))
ggplot(policies_des,aes(x=reorder(country,redis_name),y=redis_name)) +
  geom_bar(stat="identity") +
  ylab("Average funds to public benefits and services") +
  xlab("Country") +
  theme(axis.text=element_text(size=13),
        axis.title=element_text(size=17,face="bold")) +
  geom_hline(yintercept=18.25,linetype="dashed")

##########################
# Figure 2
##########################

dat$q27abin <- ifelse(dat$q27a>=3,1,0)
dat$q27bbin <- ifelse(dat$q27b>=3,1,0)
dat$q27cbin <- ifelse(dat$q27c>=3,1,0)
dat$occs <- mapvalues(dat$s27,from=c(1,2,3,4,5,6,7,8,9,10,99),
                      to=c(1,2,3,4,5,6,7,8,9,10,NA))
dat %>%
  tidyr::drop_na(q27abin,q27bbin,q27cbin) %>%
  dplyr::select(q27abin,q27bbin,q27cbin,occs) %>%
  tidyr::pivot_longer(cols=c(q27abin,q27bbin,q27cbin),
                      names_to="variable",
                      values_to="value") %>%
  dplyr::group_by(occs,variable) %>%
  dplyr::summarize(perc=100*mean(value)) %>%
  ggplot(aes(x=perc,y=reorder(occs,perc))) +
  geom_bar(aes(fill=variable),position="dodge", 
           stat="identity") +
  ylab("Occupational group") +
  xlab("Percent") + 
  scale_y_discrete(labels=c("6"="Skilled agricultural,\nforestry or fishery worker",
                            "10"="Other/\nprefer not to answer",
                            "9"="Elementary occupation",
                            "4"="Clerical support worker",
                            "3"="Technician or \nassociate professional",
                            "8"="Plant and machine operator\nor assembly worker",
                            "1"="Manager or senior offical",
                            "5"="Service or sales worker",
                            "7"="Craft or trade worker",
                            "2"="Professional")) +
  scale_fill_grey(start=.9,end=0,
                  name="Job loss due to",
                  labels=c("q27abin"="Automation",
                           "q27bbin"="Internet platform",
                           "q27cbin"="Lack of skills")) +
  theme(axis.title=element_text(size=20,face="bold"), 
        axis.text=element_text(size=16), 
        legend.title=element_text(size=20,face="bold"),
        legend.text=element_text(size=16))

##########################
# Figure 3
##########################

saut_des <- dat.des %>%
  filter(!is.na(saut)) %>%
  group_by(country) %>%
  summarize_at(vars(saut),
               list(name = mean))
ggplot(saut_des, aes(x=reorder(country,name),y=name))  +
  geom_bar(stat="identity") +
  ylab("Average subjective technological risk") +
  xlab("Country") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=20,face="bold"))

##########################
# Figure 4
##########################

dat.des2 <- as.data.frame(dat[c("benefits1y","benefits2y",
                                "country")])
ueb_des <- dat.des2 %>%
  filter(!is.na(benefits1y)) %>%
  group_by(country) %>%
  summarize_at(vars(benefits1y,benefits2y),
               list(name = mean))
ggplot(ueb_des,aes(x=reorder(country,benefits1y_name), 
                    y=benefits1y_name))  +
  geom_bar(stat="identity") +
  ylab("Net replacement rate in unemployment") +
  xlab("Country") +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=20,face="bold")) +
  geom_point(ueb_des,mapping= 
               aes(x=reorder(country,benefits1y_name), 
                   y=benefits2y_name),
             size=5)

##########################
# Model preparations
##########################

# distribution of response variables
fam.resp <- rep("gaussian",5)
               
# priors
prior.cauchy <- list(R=list(V=diag(5)/1,nu=1),  
                     G=list(G1=list(V=diag(5)/1,nu=1,
                                    alpha.mu=rep(0,5), 
                                    alpha.V=diag(25^2,5,5))))

##########################
# Figure 5
##########################

# reduce dataset to relevant variables
vars <- c("train","retrain","redis","ubi","subs","saut","country",
          "income","techies","age","edu")
dat1 <- as.data.frame(na.omit(dat[vars]))

# scale continuous variables
morph_vars <- c("train","retrain","redis","ubi","subs","saut",
                "income","age")
morph_vars_sc <- paste(morph_vars,"s",sep="_")
for (i in 1:length(morph_vars)) {
  dat1[[morph_vars_sc[i]]] <- c(rescale(dat1[[morph_vars[i]]]))
}

# run model
fmla.fig5  <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:saut_s + trait:techies + trait:age_s + trait:edu + trait:income_s  
mmlm.fig5 <- MCMCglmm(fmla.fig5,
                      random=~us(trait):country, 
                      rcov=~us(trait):units, 
                      prior=prior.cauchy,
                      data=dat1, 
                      family=fam.resp, 
                      nitt=20000,burnin=5000,thin=10)

# produce figure
coef.fig5 <- tidy(mmlm.fig5,effect="fixed") 
brackets <- list(c("TechRisk","traittrain_s:saut_s","traitsubs_s:saut_s"),
                 c("TechUsers","traittrain_s:techies1","traitsubs_s:techies1"),
                 c("Age","traittrain_s:age_s","traitsubs_s:age_s"),
                 c("HighEduc","traittrain_s:edu1","traitsubs_s:edu1"),
                 c("Income","traittrain_s:income_s","traitsubs_s:income_s"))
{dwplot(tt2,dot_args=list(size=3)) +
    geom_vline(xintercept=0,lty=2) +
    scale_y_discrete(labels=c("traittrain_s:saut_s"="Education",
                              "traitretrain_s:saut_s"="Retraining",
                              "traitredis_s:saut_s"="Benefits",
                              "traitubi_s:saut_s"="UBI",
                              "traitsubs_s:saut_s"="Subsidies",
                              "traittrain_s:techies1"="Education",
                              "traitretrain_s:techies1"="Retraining",
                              "traitredis_s:techies1"="Benefits",
                              "traitubi_s:techies1"="UBI",
                              "traitsubs_s:techies1"="Subsidies",
                              "traittrain_s:age_s"="Education",
                              "traitretrain_s:age_s"="Retraining",
                              "traitredis_s:age_s"="Benefits",
                              "traitubi_s:age_s"="UBI",
                              "traitsubs_s:age_s"="Subsidies",
                              "traittrain_s:edu1"="Education",
                              "traitretrain_s:edu1"="Retraining",
                              "traitredis_s:edu1"="Benefits",
                              "traitubi_s:edu1"="UBI",
                              "traitsubs_s:edu1"="Subsidies",
                              "traittrain_s:income_s"="Education",
                              "traitretrain_s:income_s"="Retraining",
                              "traitredis_s:income_s"="Benefits",
                              "traitubi_s:income_s"="UBI",
                              "traitsubs_s:income_s"="Subsidies")) +
    theme(axis.text=element_text(size=16)) +
    scale_color_grey()} %>%
  brack(brackets,face="bold")

##########################
# Table 2: Unemployment compensation after one year
##########################

# reduce dataset to relevant variables
vars2 <- c("train","retrain","redis","ubi","subs","saut","country",
           "income","techies","age","edu","benefits1y")
dat.ueb <- as.data.frame(na.omit(dat[vars2]))

# scale continuous variables
morph_vars2 <- c("train","retrain","redis","ubi","subs","saut",
                 "income","age","benefits1y")
morph_vars_sc2 <- paste(morph_vars2,"s",sep="_")
for (i in 1:length(morph_vars2)) {
  dat.ueb[[morph_vars_sc2[i]]] <- c(rescale(dat.ueb[[morph_vars2[i]]]))
}

# run model
fmla.tab2y1 <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:saut_s + trait:techies + trait:age_s + trait:edu + trait:income_s + trait:benefits1y_s
mmlm.tab2y1 <- MCMCglmm(mmlm.tab2y1,
                        random=~us(trait):country, 
                        rcov=~us(trait):units,
                        prior=prior.cauchy,
                        data=dat.ueb, 
                        family=fam.resp, 
                        nitt=20000,burnin=5000,thin=10)
summary(mmlm.tab2y1)

##########################
# Table 2: Unemployment compensation after two years
##########################

# reduce dataset to relevant variables
vars3 <- c("train","retrain","redis","ubi","subs","saut","country",
          "income","techies","age","edu","benefits2y")
dat.ueb2 <- as.data.frame(na.omit(dat[vars3]))

# scale continuous variables
morph_vars3 <- c("train","retrain","redis","ubi","subs","saut",
                 "income","age","benefits2y")
morph_vars_sc3 <- paste(morph_vars3,"s",sep="_")
for (i in 1:length(morph_vars3)) {
  dat.ueb2[[morph_vars_sc3[i]]] <- c(rescale(dat.ueb2[[morph_vars3[i]]]))
}

# run model
fmla.tab2y2 <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:saut_s + trait:saut_s + trait:techies + trait:age_s + trait:edu + trait:income_s + trait:benefits2y_s
mmlm.tab2y2 <- MCMCglmm(fmla.tab2y2,
                        random=~us(trait):country, 
                        rcov=~us(trait):units,
                        prior=prior.cauchy,
                        data= dat.ueb2, 
                        family=fam.resp, 
                        nitt=20000,burnin=5000,thin=10)
summary(mmlm.tab2y2)

##########################
# Figure 6
##########################

# run models
fmla.fig6y1 <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:(saut_s*benefits1y_s) + trait:techies + trait:age_s + trait:edu + trait:income_s 
mmlm.fig6y1 <- MCMCglmm(fmla.fig6y1,
                        random=~us(trait):country, 
                        rcov=~us(trait):units,
                        prior=prior.cauchy,
                        data=dat.ueb, 
                        family=fam.resp, 
                        nitt=20000,burnin=5000,thin=10)

fmla.fig6y2 <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:(saut_s*benefits2y_s) + trait:techies + trait:age_s + trait:edu + trait:income_s 
mmlm.fig6y2 <- MCMCglmm(fmla.fig6y2,
                        random=~us(trait):country, 
                        rcov=~us(trait):units,
                        prior=prior.cauchy,
                        data=dat.ueb2, 
                        family=fam.resp, 
                        nitt=20000,burnin=5000,thin=10)

# make predictions for unemployment compensation after one year
tmpdat <- as.data.frame(dat.ueb[,c("train_s","retrain_s","redis_s",
                                   "ubi_s","subs_s","saut_s","income_s",
                                   "techies","age_s","edu","country",
                                   "benefits1y_s")])
jvalues <- with(dat.ueb,seq(from=min(saut_s),to=max(saut_s),
                             length.out=10))
out=list()
for (j in jvalues[1:10]) {
  tmp_dat <- rbind(tmpdat,tmpdat)
  tmp_dat$benefits1y_s <- rep(c(-0.97121,0.75391),each=nrow(tmpdat))
  tmp_dat$saut_s <- j
  out[[length(out)+1]] <- tmp_dat
}
tmp_dat <- bind_rows(out)
pred_dat <- predict(mmlm.fig6y1,newdata=tmp_dat)

# select data for "public benefits and services" (after one year)
redis1l <- pred_dat[749561:768299]
redis1h <- pred_dat[768300:787038]
redis2l <- pred_dat[787039:805777]
redis2h <- pred_dat[805778:824516]
redis3l <- pred_dat[824517:843255]
redis3h <- pred_dat[843256:861994]
redis4l <- pred_dat[861995:880733]
redis4h <- pred_dat[880734:899472]
redis5l <- pred_dat[899473:918211]
redis5h <- pred_dat[918212:936950]
redis6l <- pred_dat[936951:955689]
redis6h <- pred_dat[955690:974428]
redis7l <- pred_dat[974429:993167]
redis7h <- pred_dat[993168:1011906]
redis8l <- pred_dat[1011907:1030645]
redis8h <- pred_dat[1030646:1049384]
redis9l <- pred_dat[1049385:1068123]
redis9h <- pred_dat[1068124:1086862]
redis10l <- pred_dat[1086863:1105601]
redis10h <- pred_dat[1105602:1124340]

redisl <- rbind(summary(redis1l),summary(redis2l),summary(redis3l),
                summary(redis4l),summary(redis5l),summary(redis6l),
                summary(redis7l),summary(redis8l),summary(redis9l),
                summary(redis10l))
redish <- rbind(summary(redis1h),summary(redis2h),summary(redis3h),
                summary(redis4h),summary(redis5h),summary(redis6h),
                summary(redis7h),summary(redis8h),summary(redis9h),
                summary(redis10h))
plotdat1.1 <- as.data.frame(cbind(redisl,jvalues))
plotdat1.2 <- as.data.frame(cbind(redish,jvalues))

# select data for "retraining" (after one year)
retr1l <- pred_dat[374781:393518]
retr1h <- pred_dat[393519:412257]
retr2l <- pred_dat[412258:430996]
retr2h <- pred_dat[430997:449735] 
retr3l <- pred_dat[449736:468474]
retr3h <- pred_dat[468475:487213]
retr4l <- pred_dat[487214:505952]
retr4h <- pred_dat[505953:524691]
retr5l <- pred_dat[524692:543430]
retr5h <- pred_dat[543431:562169]
retr6l <- pred_dat[562170:580908]
retr6h <- pred_dat[580909:599647]
retr7l <- pred_dat[599648:618386]
retr7h <- pred_dat[618387:637125]
retr8l <- pred_dat[637126:655864]
retr8h <- pred_dat[655865:674603]
retr9l <- pred_dat[674604:693342]
retr9h <- pred_dat[693343:712081]
retr10l <- pred_dat[712082:730820]
retr10h <- pred_dat[730821:749559]

retrl <- rbind(summary(retr1l),summary(retr2l),summary(retr3l),
               summary(retr4l),summary(retr5l),summary(retr6l),
               summary(retr7l),summary(retr8l),summary(retr9l),
               summary(retr10l))
retrh <- rbind(summary(retr1h),summary(retr2h),summary(retr3h),
               summary(retr4h),summary(retr5h),summary(retr6h),
               summary(retr7h),summary(retr8h),summary(retr9h),
               summary(retr10h))
plotdat1.3 <- as.data.frame(cbind(retrl, jvalues))
plotdat1.4 <- as.data.frame(cbind(retrh, jvalues))

# select data for "UBI" (after one year)
ubi1l <- pred_dat[1124341:1143079]
ubi1h <- pred_dat[1143080:1161818]
ubi2l <- pred_dat[1161819:1180557]
ubi2h <- pred_dat[1180558:1199296] 
ubi3l <- pred_dat[1199297:1218035]
ubi3h <- pred_dat[1218036:1236774]
ubi4l <- pred_dat[1236775:1255513]
ubi4h <- pred_dat[1255514:1274252]
ubi5l <- pred_dat[1274253:1292991]
ubi5h <- pred_dat[1292992:1311730]
ubi6l <- pred_dat[1311731:1330469]
ubi6h <- pred_dat[1330470:1349208]
ubi7l <- pred_dat[1349209:1367947]
ubi7h <- pred_dat[1367948:1386686]
ubi8l <- pred_dat[1386687:1405425]
ubi8h <- pred_dat[1405426:1424164]
ubi9l <- pred_dat[1424165:1442903]
ubi9h <- pred_dat[1442904:1461642]
ubi10l <- pred_dat[1461643:1480381]
ubi10h <- pred_dat[1480382:1499120]

ubil <- rbind(summary(ubi1l),summary(ubi2l),summary(ubi3l),
              summary(ubi4l),summary(ubi5l),summary(ubi6l),
              summary(ubi7l),summary(ubi8l),summary(ubi9l),
              summary(ubi10l))
ubih <- rbind(summary(ubi1h),summary(ubi2h),summary(ubi3h),
              summary(ubi4h),summary(ubi5h),summary(ubi6h),
              summary(ubi7h),summary(ubi8h),summary(ubi9h),
              summary(ubi10h))
plotdat1.5 <- as.data.frame(cbind(ubil,jvalues))
plotdat1.6 <- as.data.frame(cbind(ubih,jvalues))

# make predictions for unemployment compensation after two years
tmpdat2 <- as.data.frame(dat.ueb2[,c("train_s","retrain_s","redis_s",
                                     "ubi_s","subs_s","saut_s","income_s",
                                     "techies","age_s","edu","country",
                                     "benefits2y_s")])
jvalues2 <- with(dat.ueb2,seq(from=min(saut_s),to=max(saut_s),
                              length.out=10))
out2=list()
for (j in jvalues2[1:10]) {
  tmp_dat2 <- rbind(tmpdat2,tmpdat2)
  tmp_dat2$benefits2y_s <- rep(c(-0.87464,0.92258),each=nrow(tmpdat2))
  tmp_dat2$saut_s <- j
  out2[[length(out2) + 1]] <- tmp_dat2
}
tmp_dat2 <- bind_rows(out2)
pred_dat2 <- predict(mmlm.fig6y2,newdata=tmp_dat2)

# select data for "public benefits and services" (after two years)
redis1l <- pred_dat2[749561:768298]
redis1h <- pred_dat2[768299:787037]
redis2l <- pred_dat2[787038:805776]
redis2h <- pred_dat2[805777:824515] 
redis3l <- pred_dat2[824516:843254]
redis3h <- pred_dat2[843255:861993]
redis4l <- pred_dat2[861994:880732]
redis4h <- pred_dat2[880733:899471]
redis5l <- pred_dat2[899472:918210]
redis5h <- pred_dat2[918211:936949]
redis6l <- pred_dat2[936950:955688]
redis6h <- pred_dat2[955689:974427]
redis7l <- pred_dat2[974428:993166]
redis7h <- pred_dat2[993167:1011905]
redis8l <- pred_dat2[1011906:1030644]
redis8h <- pred_dat2[1030645:1049383]
redis9l <- pred_dat2[1049384:1068122]
redis9h <- pred_dat2[1068123:1086861]
redis10l <- pred_dat2[1086862:1105600]
redis10h <- pred_dat2[1105601:1124339]

redisl <- rbind(summary(redis1l),summary(redis2l),summary(redis3l),
                summary(redis4l),summary(redis5l),summary(redis6l),
                summary(redis7l),summary(redis8l),summary(redis9l),
                summary(redis10l))
redish <- rbind(summary(redis1h),summary(redis2h),summary(redis3h),
                summary(redis4h),summary(redis5h),summary(redis6h),
                summary(redis7h),summary(redis8h),summary(redis9h),
                summary(redis10h))
plotdat2.1 <- as.data.frame(cbind(redisl, jvalues2))
plotdat2.2 <- as.data.frame(cbind(redish, jvalues2))

# select data for "retraining" (after two years)
retr1l <- pred_dat2[374781:393518]
retr1h <- pred_dat2[393519:412257]
retr2l <- pred_dat2[412258:430996]
retr2h <- pred_dat2[430997:449735] 
retr3l <- pred_dat2[449736:468474]
retr3h <- pred_dat2[468475:487213]
retr4l <- pred_dat2[487214:505952]
retr4h <- pred_dat2[505953:524691]
retr5l <- pred_dat2[524692:543430]
retr5h <- pred_dat2[543431:562169]
retr6l <- pred_dat2[562170:580908]
retr6h <- pred_dat2[580909:599647]
retr7l <- pred_dat2[599648:618386]
retr7h <- pred_dat2[618387:637125]
retr8l <- pred_dat2[637126:655864]
retr8h <- pred_dat2[655865:674603]
retr9l <- pred_dat2[674604:693342]
retr9h <- pred_dat2[693343:712081]
retr10l <- pred_dat2[712082:730820]
retr10h <- pred_dat2[730821:749559]

retrl <- rbind(summary(retr1l),summary(retr2l),summary(retr3l),
               summary(retr4l),summary(retr5l),summary(retr6l),
               summary(retr7l),summary(retr8l),summary(retr9l),
               summary(retr10l))
retrh <- rbind(summary(retr1h),summary(retr2h),summary(retr3h),
               summary(retr4h),summary(retr5h),summary(retr6h),
               summary(retr7h),summary(retr8h),summary(retr9h),
               summary(retr10h))
plotdat2.3 <- as.data.frame(cbind(retrl,jvalues2))
plotdat2.4 <- as.data.frame(cbind(retrh,jvalues2))

# select data for "UBI" (after two years)
ubi1l <- pred_dat2[1124341:1143079]
ubi1h <- pred_dat2[1143080:1161818]
ubi2l <- pred_dat2[1161819:1180557]
ubi2h <- pred_dat2[1180558:1199296] 
ubi3l <- pred_dat2[1199297:1218035]
ubi3h <- pred_dat2[1218036:1236774]
ubi4l <- pred_dat2[1236775:1255513]
ubi4h <- pred_dat2[1255514:1274252]
ubi5l <- pred_dat2[1274253:1292991]
ubi5h <- pred_dat2[1292992:1311730]
ubi6l <- pred_dat2[1311731:1330469]
ubi6h <- pred_dat2[1330470:1349208]
ubi7l <- pred_dat2[1349209:1367947]
ubi7h <- pred_dat2[1367948:1386686]
ubi8l <- pred_dat2[1386687:1405425]
ubi8h <- pred_dat2[1405426:1424164]
ubi9l <- pred_dat2[1424165:1442903]
ubi9h <- pred_dat2[1442904:1461642]
ubi10l <- pred_dat2[1461643:1480381]
ubi10h <- pred_dat2[1480382:1499120]

ubil <- rbind(summary(ubi1l),summary(ubi2l),summary(ubi3l),
              summary(ubi4l),summary(ubi5l),summary(ubi6l),
              summary(ubi7l),summary(ubi8l),summary(ubi9l),
              summary(ubi10l))
ubih <- rbind(summary(ubi1h),summary(ubi2h),summary(ubi3h),
              summary(ubi4h),summary(ubi5h),summary(ubi6h),
              summary(ubi7h),summary(ubi8h),summary(ubi9h),
              summary(ubi10h))
plotdat2.5 <- as.data.frame(cbind(ubil,jvalues2))
plotdat2.6 <- as.data.frame(cbind(ubih,jvalues2))

# produce graph
par(mfrow=c(3,2),cex=1.2)
plot(plotdat1.3$jvalues, plotdat1.3$Mean,type ="l",lwd=2,
     cex=2,
     ylim=c(-0.07,0.17), 
     ylab="Funding of retraining",
     xlab="",
     main="(a) Unemployment compensation after 1 year: \nRetraining",
     panel.first=c(abline(h=seq(-0.05,0.15,0.05),lty=1,col="#d9d9d9"),
                   abline(v=seq(-0.5,0.5,0.5),lty=1,col="#d9d9d9")))
polygon(c(plotdat1.3$jvalues,rev(plotdat1.3$jvalues)),
        c(plotdat1.3$`1st Qu.`,rev(plotdat1.3$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f = 0.3),border=NA)
polygon(c(plotdat1.4$jvalues,rev(plotdat1.4$jvalues)),
        c(plotdat1.4$`1st Qu.`,rev(plotdat1.4$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f=0.3),border=NA)
lines(plotdat1.4$jvalues, plotdat1.4$Mean, lwd = 2, lty = 2)
legend("topright",c("Low compensation","High compensation"),
       lty=c(1,2),bty="n",lwd=c(2,2))
plot(plotdat2.3$jvalues2,plotdat2.3$Mean,type="l",lwd = 2,
     ylim=c(-0.07,0.17), 
     ylab="",
     xlab="",
     main="(b) Unemployment compensation after 2 years: \nRetraining",
     panel.first=c(abline(h=seq(-0.05,0.20,0.05),lty=1,col="#d9d9d9"),
                   abline(v=seq(-0.5,0.5,0.5),lty=1,col="#d9d9d9")))
polygon(c(plotdat2.3$jvalues,rev(plotdat2.3$jvalues2)),
        c(plotdat2.3$`1st Qu.`,rev(plotdat2.3$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f=0.3),border=NA)
polygon(c(plotdat2.4$jvalues,rev(plotdat2.4$jvalues2)),
        c(plotdat2.4$`1st Qu.`,rev(plotdat2.4$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f=0.3),border= NA)
lines(plotdat2.4$jvalues2,plotdat2.4$Mean,lwd = 2,lty=2)
legend("topright",c("Low compensation","High compensation"),
       lty=c(1,2),bty ="n",lwd =c(2,2))
plot(plotdat1.1$jvalues,plotdat1.1$Mean,type="l",lwd=2,
     ylim=c(-0.2,0.1),yaxt ="n",
     ylab="Funding of public benefits and services",
     xlab="",
     main="Public benefits and services",
     panel.first=c(abline(h=seq(-0.20,0.10,0.10),lty=1,col="#d9d9d9"),
                   abline(v=seq(-0.5,0.5,0.5),lty=1,col="#d9d9d9")))
axis(2,at=seq(-0.2,0.1,by=0.1))
polygon(c(plotdat1.1$jvalues,rev(plotdat1.1$jvalues)),
        c(plotdat1.1$`1st Qu.`,rev(plotdat1.1$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f = 0.3),border=NA)
polygon(c(plotdat1.2$jvalues,rev(plotdat1.2$jvalues)),
        c(plotdat1.2$`1st Qu.`,rev(plotdat1.2$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f=0.3),border=NA)
lines(plotdat1.2$jvalues,plotdat1.2$Mean,lwd=2,lty=2)
legend("bottomright",c("Low compensation","High compensation"),
       lty=c(1,2),bty ="n",lwd=c(2,2))
plot(plotdat2.1$jvalues2,plotdat2.1$Mean,type = "l",lwd=2,
     ylim=c(-0.2,0.1),yaxt="n",
     ylab="",
     xlab="",
     main="Public benefits and services",
     panel.first=c(abline(h=seq(-0.20,0.10,0.10),lty=1,col="#d9d9d9"),
                   abline(v=seq(-0.5,0.5,0.5),lty=1,col="#d9d9d9")))
axis(2,at=seq(-0.20,0.10,by=0.10))
polygon(c(plotdat2.1$jvalues2,rev(plotdat2.1$jvalues2)),
        c(plotdat2.1$`1st Qu.`,rev(plotdat2.1$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f=0.3), border=NA)
polygon(c(plotdat2.2$jvalues,rev(plotdat2.2$jvalues2)),
        c(plotdat2.2$`1st Qu.`,rev(plotdat2.2$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f = 0.3),border=NA)
lines(plotdat2.2$jvalues2,plotdat2.2$Mean,lwd=2,lty=2)
legend("bottomright",c("Low compensation","High compensation"),
       lty=c(1,2),bty ="n",lwd=c(2,2))
plot(plotdat1.5$jvalues,plotdat1.5$Mean,type="l",lwd=2,
     ylim=c(-0.10,0.05), 
     ylab="Funding of UBI",
     xlab="Subjective technological risk",
     main="UBI",
     panel.first=c(abline(h=seq(-0.05,0.20,0.05),lty=1,col="#d9d9d9"),
                   abline(v=seq(-0.5,0.5,0.5),lty=1,col="#d9d9d9")))
polygon(c(plotdat1.5$jvalues,rev(plotdat1.5$jvalues)),
        c(plotdat1.5$`1st Qu.`,rev(plotdat1.5$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f = 0.3),border=NA)
polygon(c(plotdat1.6$jvalues,rev(plotdat1.6$jvalues)),
        c(plotdat1.6$`1st Qu.`,rev(plotdat1.6$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f=0.3),border=NA)
lines(plotdat1.6$jvalues,plotdat1.6$Mean,lwd=2,lty=2)
legend("bottomright",c("Low compensation","High compensation"),
       lty=c(1,2),bty="n",lwd=c(2,2))
plot(plotdat2.5$jvalues2,plotdat2.5$Mean,type="l",lwd=2,
     ylim=c(-0.10,0.05), 
     ylab="",
     xlab="Subjective technological risk",
     main="UBI",
     panel.first=c(abline(h=seq(-0.05,0.20,0.05),lty=1,col="#d9d9d9"),
                   abline(v=seq(-0.5,0.5,0.5),lty=1,col="#d9d9d9")))
polygon(c(plotdat2.5$jvalues2,rev(plotdat2.5$jvalues2)),
        c(plotdat2.5$`1st Qu.`,rev(plotdat2.5$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f = 0.3),border=NA)
polygon(c(plotdat2.6$jvalues2,rev(plotdat2.6$jvalues2)),
        c(plotdat2.6$`1st Qu.`,rev(plotdat2.6$`3rd Qu.`)),
        col=adjustcolor("grey",alpha.f = 0.3),border=NA)
lines(plotdat2.6$jvalues2, plotdat2.6$Mean, lwd = 2, lty = 2)
legend("bottomright",c("Low compensation","High compensation"),
       lty=c(1,2),bty="n",lwd=c(2,2))

##########################
# Figure 7
##########################

# RTI scores
# reduce dataset to relevant variables
varsrti <- c("train","retrain","redis","ubi","subs","rti_score",
             "country","income","techies","age","edu")
dat.rti <- as.data.frame(na.omit(dat[varsrti]))
# scale continuous variables
morph_varsrti <- c("train","retrain","redis","ubi","subs","rti_score",
                   "income","age")
morph_vars_scrti <- paste(morph_varsrti,"s",sep="_")
for (i in 1:length(morph_varsrti)) {
  dat.rti[[morph_vars_scrti[i]]] <- c(rescale(dat.rti[[morph_varsrti[i]]]))
}
# run model 
fmla.fig7rti  <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:rti_score_s + trait:techies + trait:age_s + trait:edu + trait:income_s 
mmlm.fig7rti <- MCMCglmm(fmla.fig7rti,
                         random=~us(trait):country, 
                         rcov=~us(trait):units,
                         prior=prior.cauchy,
                         data=dat.rti, 
                         family=fam.resp, 
                         nitt=20000,burnin=5000,thin=10)
coef.fig7rti <- tidy(mmlm.fig7rti,effect="fixed")

# individual components of subjective technological risk
# reduce dataset to relevant variables
varsq27 <- c("train","retrain","redis","ubi","subs","q27a","q27b","q27c",
             "country","income","techies","age","edu")
dat.q27 <- as.data.frame(na.omit(dat[varsq27]))

# scale continuous variables
morph_varsq27 <- c("train","retrain","redis","ubi","subs","q27a","q27b",
                   "q27c","income","age")
morph_vars_scq27 <- paste(morph_varsq27,"s",sep="_")
for (i in 1:length(morph_varsq27)) {
  dat.q27[[morph_vars_scq27[i]]] <- c(rescale(dat.q27[[morph_varsq27[i]]]))
}
# first component
fmla.fig7q27a  <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:q27a_s + trait:techies + trait:age_s + trait:edu + trait:income_s  
mmlm.fig7q27a <- MCMCglmm(fmla.fig7q27a,
                          random=~us(trait):country, 
                          rcov=~us(trait):units,
                          prior= prior.cauchy,
                          data=dat.q27, 
                          family=fam.resp, 
                          nitt=20000,burnin=5000,thin=10) 
coef.fig7q27a <- tidy(mmlm.fig7q27a,effect ="fixed")
# second component
fmla.fig7q27b  <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:q27b_s + trait:techies + trait:age_s + trait:edu + trait:income_s  
mmlm.fig7q27b <- MCMCglmm(fmla.fig7q27b,
                          random=~us(trait):country, 
                          rcov=~us(trait):units,
                          prior=prior.cauchy,
                          data=dat.q27,
                          family=fam.resp, 
                          nitt=20000,burnin=5000,thin=10) 
coef.fig7q27b <- tidy(mmlm.fig7q27b,effect="fixed")
# third component
fmla.fig7q27c  <- cbind(train_s,retrain_s,redis_s,ubi_s,subs_s) ~
  -1 + trait:q27c_s + trait:techies + trait:age_s + trait:edu + trait:income_s  
mmlm.fig7q27c <- MCMCglmm(fmla.fig7q27c,
                         random=~us(trait):country, 
                         rcov=~us(trait):units,
                         prior=prior.cauchy,
                         data=dat.q27, 
                         family=fam.resp, 
                         nitt=20000,burnin=5000,thin=10)
coef.fig7q27c <- tidy(mmlm.fig7q27c,effect="fixed")

# produce graph
coef.fig7 <- bind_rows(RTI=coef.fig7rti,Subjective.3=coef.fig7q27c, 
                       Subjective.2=coef.fig7q27b, 
                       Subjective.1=coef.fig7q27a, 
                       Subjective=coef.fig5,.id="model") %>%
  filter(effect=="fixed")
coef.fig7$term[1] <- c("traittrain_s:saut_s")
coef.fig7$term[2] <- c("traitretrain_s:saut_s")
coef.fig7$term[3] <- c("traitredis_s:saut_s")
coef.fig7$term[4] <- c("traitubi_s:saut_s")
coef.fig7$term[5] <- c("traitsubs_s:saut_s")
coef.fig7$term[26] <- c("traittrain_s:saut_s")
coef.fig7$term[27] <- c("traitretrain_s:saut_s") 
coef.fig7$term[28] <- c("traitredis_s:saut_s")
coef.fig7$term[29] <- c("traitubi_s:saut_s")
coef.fig7$term[30] <- c("traitsubs_s:saut_s")
coef.fig7$term[51] <- c("traittrain_s:saut_s")
coef.fig7$term[52] <- c("traitretrain_s:saut_s")
coef.fig7$term[53] <- c("traitredis_s:saut_s")
coef.fig7$term[54] <- c("traitubi_s:saut_s")
coef.fig7$term[55] <- c("traitsubs_s:saut_s")
coef.fig7$term[76] <- c("traittrain_s:saut_s")
coef.fig7$term[77] <- c("traitretrain_s:saut_s")
coef.fig7$term[78] <- c("traitredis_s:saut_s")
coef.fig7$term[79] <- c("traitubi_s:saut_s")
coef.fig7$term[80] <- c("traitsubs_s:saut_s")

{dwplot(coef.fig7,dodge_size=0.8,dot_args=list(aes(shape = model),
                                               size = 2.5)) +
    geom_vline(xintercept=0,lty=2) +
    scale_y_discrete(labels=c("traittrain_s:saut_s"="Education",
                              "traitretrain_s:saut_s"="Retraining",
                              "traitredis_s:saut_s"="Benefits",
                              "traitubi_s:saut_s"="UBI",
                              "traitsubs_s:saut_s"="Subsidies",
                              "traittrain_s:techies1"="Education",
                              "traitretrain_s:techies1"="Retraining",
                              "traitredis_s:techies1"="Benefits",
                              "traitubi_s:techies1"="UBI",
                              "traitsubs_s:techies1"="Subsidies",
                              "traittrain_s:age_s"="Education",
                              "traitretrain_s:age_s"="Retraining",
                              "traitredis_s:age_s"="Benefits",
                              "traitubi_s:age_s"="UBI",
                              "traitsubs_s:age_s"="Subsidies",
                              "traittrain_s:edu1"="Education",
                              "traitretrain_s:edu1"="Retraining",
                              "traitredis_s:edu1"="Benefits",
                              "traitubi_s:edu1"="UBI",
                              "traitsubs_s:edu1"="Subsidies",
                              "traittrain_s:income_s"="Education",
                              "traitretrain_s:income_s"="Retraining",
                              "traitredis_s:income_s"="Benefits",
                              "traitubi_s:income_s"="UBI",
                              "traitsubs_s:income_s"="Subsidies")) +
    theme(axis.text=element_text(size=16), 
          legend.title=element_text(size=16),
          legend.text=element_text(size=16),
          legend.position=c(1, 0.907),
          legend.background=element_rect(colour="grey80")) +
    scale_color_grey(name="Measure of TechRisk",
                     breaks=c("Subjective","Subjective.1","Subjective.2",
                              "Subjective.3","RTI")) +
    scale_shape_discrete(name="Measure of TechRisk",
                         breaks=c("Subjective","Subjective.1",
                                  "Subjective.2","Subjective.3",
                                  "RTI"))} %>%
  brack(brackets,face="bold")  
   